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We examine the star lattice Kitaev model whose ground state is a a chiral spin liquid. We fermionize the model 
such that the fermionic vacua are toric code states on an effective Kagome lattice. This implies that the Abelian 
phase of the system is inherited from the fermionic vacua and that time reversal symmetry is spontaneously 
broken at the level of the vacuum. In terms of these fermions we derive the Bloch-matrix Hamiltonians for the 
vortex free sector and its time reversed counterpart and illuminate the relationships between the sectors. The 
phase diagram for the model is shown to be a sphere in the space of coupling parameters around the triangles 
of the lattices. The Abelian phase lies inside the sphere and the critical boundary between topologically distinct 
Abelian and non- Abelian phases lies on the surface. Outside the sphere the system is generically gapped except 
in the planes where the coupling parameters between the vertices on triangles are zero. These cases correspond 
to bipartite lattice structures and the dispersion relations are similar to that of the original Kitaev honeycomb 
model. In a further analysis we demonstrate the three-fold non- Abelian groundstate degeneracy on a torus by 
explicit calculation. 
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I. INTRODUCTION 

In his original analysis of the honeycomb model 1 , Kitaev 
noted that a similar type of system but with triangles placed 
at the vertices of the honeycomb lattice would spontaneously 
break time reversal symmetry. A system of precisely this type 
was subsequently analyzed by Yao and Kivelson 2 and shown 
to be an example of a chiral spin liquid with just nearest neigh- 
bor interactions between sites. This system also inherits a 
number of interesting properties from the original honeycomb 
model, particularly the existence of both Abelian and non- 
Abelian topological phases. A finite temperature analysis of 
the model has recently been performed 3 and, as with the orig- 
inal honeycomb model, there are important overlaps with the 
physics of classical dimer models and Kasteleyn matrices^. 

In what follows we present an analysis of this system using 
the Jordan- Wigner fermionization procedure of Ref. 5 which 
explicitly formulates the fermionic vacua as toric code states 
on an effective Kagome lattice. The fermionization procedure 
is a two step process where we map the model to a system of 
hard-core bosons and spins on an effective Kagome lattice^El 
and then define fermions in terms of the hard-core bosons and 
spins. Once fermionized, and similarly to the original honey- 
comb lattice model 8 , the ground-state sector of the system can 
be transformed to that of a spinless p-wave superconductor. 
With our method however, we obtain vacua for the fermion- 
ized theory as the stabilized wavefunctions of an abelian toric 
code model 9 , defined with effective spins on a Kagome lat- 
tice. These vacuum wavefunctions are independent of the cou- 
plings of the model. The ground-state for the full system, valid 
for the entire parameter space, is a BCS type condensate over 
the toric code ground-state. The topological degeneracies of 
the model are already present at the level of the fermionic vac- 
uum. However, by generalizing the arguments presented irJ5U2l 
we show how the BCS product lifts some of this degeneracy in 
the non- Abelian phase. The predicted degeneracy is in agree- 
ment with the original observations of Yao and KivelsorP. 



The representation we use also illuminates the nature of the 
spontaneously broken time-reversal symmetry. We see for ex- 
ample that this symmetry is broken at the level of the vacua. 
The chiral nature of these vacuum states has been recently 
analyzed 11 . A detailed analysis of the phase boundary be- 
tween the Abelian and the non- Abelian phase is also included 
and we show that this critical boundary can be understood as 
a sphere in the space of the spin-spin coupling strengths. We 
also see that the opening of a gap to a non-Abelian phase is 
highly dependent on the interaction strengths around the trian- 
gular plaquettes. We will show that if any of these interactions 
are subtracted from the Hamiltonian we cannot open a gap to 
the non-abelian phase. Exact expressions for the energy dis- 
persion relations for ground states sectors are proved in the 
appendix, which also provides a brief analysis of the Dirac 
cone structures of the gapless phases. 



II. THE STAR LATTICE KITAEV MODEL AND THE 
HEXAGONAL TORIC CODE 

The Hamiltonian consists of directional spin-spin interac- 
tions on the star lattice (also known as the triangle-honeycomb 
or Fisher lattice). We use the representation of the model in- 
troduced in Ref. 6 as it provides a straightforward route , by 
contracting the Z-links, to the definition of the fermions as 
toric code states on the Kagome lattice. In this representation 
the Hamiltonian can be written as 



H 



H z + Hj + H K + H L 

-z aZ ° z - J ^~ 



(1) 



Z-links 



J-links 



x rrV 



K 

K-links 



x rrV 



L °~n 

L-links 



where it should be understood that the Z-links connect sepa- 
rate triangles and the J,K and L-links within the triangles are 
the positive, zero and negative slopes respectively, see Figure 
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FIG. 1: (Color online) The star lattice is a hexagonal lattice with 
each vertex replaced by a triangle. There are two types of plaque- 
ttes symmetries in the system: triangular and dodecagonal. The tri- 
angular symmetries are responsible for the spontaneous breaking of 
time-reversal symmetry-^ 




FIG. 2: (Color online) The effective Kagome lattice and plaquette 
operators. Note that within our labeling convention the A q operator 
does not actually contain any sites from the unit cell at q. 



[T] We refer to triangles that point up as black triangles and 
those that point down as white triangles. The sites on these 
triangles in the original lattice are colored black and white re- 
spectively. 

We define a basic unit cell of the lattice around a white tri- 
angle. We label the Z-link at the bottom of the triangle with 
n = 1, the Z-link from the top right with n = 2 and the Z-link 
from the top left with n = 3. Each spin site can be specified 
using the position vector q, the index n, and whether it is on a 
□ or a ■ site. In a 6iV spin system we have N unit cells. 

Using this basic unit cell we we may write 
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where have introduced here the shorthand notation 
q \= q — n x + n y , q |= q — n y etc. and the two unit 
vectors n x and n y as shown in FIG. [T] We have set all 
coupling strengths on all z-dimers to Z thus restricting the 
parameter space of the original model in this direction. 

Within the model, as in honeycomb lattice model, there are 
closed loop symmetries that we can generate (up to a phase) 
with overlapping link operators. The simplest of these are the 
dodecagonal and triangular loops. These are defined pictori- 
ally in Figure [T] For simplicity we will refer to a generic loop 
symmetry as W q and define it such that it has eigenvalues of 
±1. The fact that the Hamiltonian commutes with all plaque- 
tte operators implies that we may choose energy eigenvectors 



n) such that W q = (n \W q \ 



= ±1. lfW q = -1 then we 



say that the state | n) carries either a triangular or dodecago- 
nal vortex at q. When we refer to a vortex-sector we mean the 
subspace of the system Hilbert space with a particular con- 
figuration of vortices. The vortex-free sector for example is 
the subspace spanned by all eigenvectors such that W q = 1 
for all q. We will have cause later to distinguish between tri- 
angular and dodecagonal plaquettes. Our labeling convention 
will however reflect the Kagome lattice on which the fermions 
(and hard-core bosons and effective spins) are defined. 

On a torus, the plaquette operators are not independent, 
as they obey Y[ W q = I where the product is over all 
q. There are also two independent homologically non-trivial 
loop symmetries. To represent these we are free to choose 
any two closed loop operators that traverse the torus as 
long as they cannot be deformed into each other by plaque- 
tte multiplication. All other homologically non-trivial loop 
symmetries can be constructed from the products of these 
two operators and the 37V — 1 independent plaquette op- 
erators. When the torus is specified by periodic boundary 
vectors (x, y) which are integer multiples of the unit vec- 
tors i.e. x — N x n x and y = N y n y , it is natural to 
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non-trivial symmetries. We will generally use the operators 



L { x) and L y Q y ' that run through the origin as the two indepen- 
dent symmetries. For an analysis of the loop symmetries in 
the original honeycomb model see Ref . [12] 

It was pointed out by Kitaev that a model of this type (i.e 
with triangles at the vertices of a honeycomb lattice) must 
spontaneously break time reversal symmetry. This implies 
that all states of the system are at least two fold degenerate 1 . 
To see this let T be the time reversal operator. Now, because 
Ta a = — er Q , the time reversal operator changes the eigenval- 
ues of all triangular plaquette operators. However, the opera- 
tor itself commutes with the Hamiltonian as it contains only 
terms of the form a a a b . Each eigenstate must therefore have a 
time-reversed counterpart with the same energy, but from the 
vortex sector where the eigenvalues of all triangular plaquettes 
are negated^. 

The Hamiltonian ([T| can be written in terms of hard- 
core bosons and effective spins of the Z-dimers using the 
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The labels on the left hand side indicate the states of the 
z-dimer in the S z basis. The first quantum number of the 
kets on the right hand side represents the effective spin of the 
square lattice and the second is the hard-core bosonic occu- 
pation number. The presence of a boson indicates an anti- 
ferromagnetic configuration of the spins connected by a z- 
link. 

In the Abelian phase, the dominance of the Z-coupling 
terms means that spins on these links tend to align. In this 
limit the bosons are energetically suppressed. A perturba- 
tive analysis for the low energy effective Hamiltonian in this 
regime shows that the first non-constant terms occur at the 6th 
order and 8th orders. We have from Ref. 6 

ffeff = -Jg Y, °1~ J eff E P ^O q , V qfc ) (7) 

q q 

With Oq = V, = ftp T q and *q = II. T q Where 

7"° is the Pauli operator acting on the effective spin at posi- 
tion q. The functional P refers to combinations of a hexagon 
with two attached triangles. This effective Hamiltonian, de- 
fined now on a Kagome lattice, is unitarily equivalent to what 
is known as the hexagonal toric code, see for example Ref-fTTl 
The hexagonal toric code shares many of the same properties 
as the original square toric code system of Kitaev 9 . All eigen- 
states of the effective Hamiltonian (|7]i on a plane may be com- 
pletely characterized by the set of eigenvalues {O q , V q , k q } 
for all q. The ground states are those with all plaquette eigen- 
values equal to +1 and its time reversed counterpart with all 
triangular plaquettes equal to -1 and on all hexagonal plaque- 
tte eigenvalues equal to +1 2 6 . On a torus of N unit cells with 
3N effective spins we have the following identities: 



N 



N 



I. 



I 



(8) 



and so we have there a total of 3N — 2 independent plaquette 
symmetries. However we gain two independent homologi- 
cally non-trivial symmetries and thus eigenstates on a torus 
are uniquely labeled by using the full set of independent sym- 
metries. 

The basis (|6| also describes anti-ferromagnetic configura- 
tions of the z-dimers through the bosonic occupation number 
and forms an orthonormal basis for the full star lattice system. 
The Pauli operators of the original spin Hamiltonian can be 
written as (see Refs.|7]l : 



= r^{b\ + b q ) , al a =b\ + b q , 

= TV(bl + b q ) , ay q<n =iT*(b q -b q ), 
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where and b are the creation and annihilation operators for 
the hard-core bosons. In this representation we have 



(10) 




FIG. 3: (Color online) The three strings S Qt i, S Ql 2, S q ,3 in the effec- 
tive spin hard-core boson representation. Here a — (2b^b — I). 
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for the triangles. The unit cell and the plaquettes are shown 
pictorially in Figure|2] Note that with this labeling convention 
the k q operator does not contain any sites from within the unit 
cell at q. 



III. FERMIONIZATION 

To fermionize the problem we follow the Jordan- Wigner 
procedure used in Ref. [5] for the honeycomb lattice model. 
To each of the hard-core bosons in the unit cell we attach a 
string. The strings for hard-core bosons inside the same unit 
cell are identical everywhere except inside the unit cell, where 
they branch out. These strings, shown pictorially in Fig [3] are 
designed so that they break/fix a Z-dimer at the desired point 
on the lattice. This will ensure that the presence of a fermion 
implies an anti-ferromagnetic configuration of the z-dimer. 
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FIG. 4: (Color online)The operator D q and the vector q. The oper- 
ator D q is the product of all hexagonal and triangular plaquettes in 
the shaded region 



Note that the start of each string S q violates the plaquette 
symmetries Aq and Oo^- This linking of fermions and vortex- 
pairs seems to be a generic property of Kitaev models that has 
important consequences for the spectral properties on a torus 5 

Substituting the 6 y n = ff' 4,n and 6 q „ = S^ n c qn 
into expressions (10) , ( [TT) , ( fl"2) and ( fl"3j ) gives the following 
fermionic expressions for the Hamiltonian terms: 

H z = zY / (KnC q , n ~I) (21) 
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We define the string operators using the following conven- 
tion: First apply a single o\ 1 term to a black site of the Z-link 
which we set to be the origin. The rest of the string is made by 
operating cyclically with the overlapping links of the Hamil- 
tonian, ofcrj, (Tj&k' a k a i an( ^ <J i <J m^ m tne n x direction until 
we reach a required length and then a^a^, <j'^o y , a z u^ and 
VpCq in the n y direction until we reach the unit cell at q. The 
ends of the string depend on the Z-link in question. 

In the original spin representation we can write 
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The hard-core-boson/effective-spin representation of the 
strings for each of the three Z-links is shown in Figure[3] Note 
that the end of the string always contains a term (b^ + b ) or 
(&t — b ) thus ensuring a hard-core boson is created or de- 
stroyed. 

The operators S q ^ n square to unity while different operators 
Sq, m ,S q i jn anti-commute with each other. This leads us to 
identify the strings S qy „ with the following sum of fermionic 
creation and annihilation operators: S qn = c q n + c q n = 

(b q .n + b q>n )S qn where S q n is simply the string S q ^ n but 
with the (b q n + b q n ) bosonic dependence of the end-point 
removed, see Figure [5] Individually our fermionic canonical 
creation and annihilation operators are 

C q,n = bq :n S q n , c q,n = b q . n S q n (19) 

where the strings now ensure that all operators c q n and c q .„ 
obey the canonical fermionic anti -commutator relations 
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where □ is the rectangular product of plaquette operators 
(hexagons and triangles) shown in the Figure |4] On a torus the 
terms that connect opposite sides will have some dependence 
on the homologically non-trivial loop symmetries. Details on 
how to calculate their precise values can be found in Ref. |5] 

The Jordan-Wigner transform has been chosen so that the 
vacuum states for the fermions are the toric code states on the 
effective Kagome lattice. This immediately implies that that 
under time reversal the fermionic vacua must be exchanged. 
The form of the fermionic Hamiltonian also reveals a number 
of important features. We see that the triangular vorticity is 
incorporated within the Hj term and the Hk term, through 
the D q operator. However, the eigenvalue of the D q does not 
change under time reversal as it contains an even number of 
triangles. This means that if K = then eigenstates | ip) 
and T\ ip) are fermionically identical. This would seem like 
a contradiction ( recall these are different sectors ) but we are 
saved by recalling that our formulation also demands that the 
vacuum sectors be defined in terms of the hexagonal toric code 
on the effective Kagome lattice. We see that fermionically 
these eigenstates are the same, they have the same fermion 
density and the same fermion number parity. Indeed they are 
structurally identical in every way except for the vacuum from 
which they were created. In the opposite sense, in terms of the 
fermions at least, any sign of spontaneous symmetry breaking 
only occurs in the terms in Hk, which always contain an odd 
number of triangular plaquettes. These terms closely resemble 
the time reversal invariant terms in the honeycomb modef^. 
We will see later that, in the same way, these terms are jointly 
responsible for the opening of a gap. 

To proceed we first re-write the fermonic Hamiltonian as 
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where the q now label both the position and internal indices. 
The system is diagonalized by solving the Bogoliubov-De 
Gennes eigenvalue problem 

t 

, (26) 

where the non-zero entries of the diagonal matrix E nm = 
E n Snm are the quasiparticle excitation energies. The 
Bogoliubov- Valentin quasiparticle excitations are 
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Inverting ( 27 1 and substituting into ( 25 i gives 
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Normally one assumes that all the values of E n in this equa- 
tion are positive. It is this choice that one usually uses to 
obtain the ground state energy of — ^ E n /2. However, the 
Hartree-Fock-Bogoliubov formulation above actually only re- 



quires that the values of E n in (25 1 come as negated pairs. 
It does not specify that positive energies must be associated 
with 7T operators rather than 7 operators. This is an important 
point as physical situations arise naturally where the annihila- 
tion of a quasi-particle costs energy. 

In the vortex-free sector (eigenvalues of all plaquettes are 
+ 1) we can move to momentum space with the Fourier Trans- 
form 
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After some manipulation we then arrive at the following mo- 
mentum space representation for the planar Hamiltonian 
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FIG. 5: (Color online) The minimum energy gap of the vortex free 
sector with Z — 1. The critical point can be observed along the 
\J J 2 + K 2 + L 2 — 1 line. The system is gapless when J — L > 
1 / v2 and K = 0. More generally if any of the parameters J, K, or 
L are zero the system is gapless beyond the phase transition. 




FIG. 6: (Color online) Schematic of the system phase diagram. The 
surface of the sphere of radius Z indicates the critical boundary be- 
tween abelian and non-abelian phases. Inside the sphere we have a 
gapped abelian phase. Outside the sphere we are in a gapped non- 
abelian phase, provided we are not on the J = 0, K = 01 L = 
planes indicated in light green. On these planes the system is gapless. 
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and we have set 9 X = exp(ik x ) and y — exp(ik y ). 

In position space the eigenvalues of the BdG equation come 
in positive and negative pairs with eigenvectors of the form 
(U, V) T and (V*,U*) T respectively. In momentum repre- 
sentation however the positive and negative eigenvectors take 
the form (U, V) T and (— PV* , PU*) T respectively where the 
3x3 diagonal matrix P has elements Pn = 1, P22 = e" 
and P 33 = e" 



From ( 23 ) we see that the effect of the time reversal oper- 
ator T is to change the signs of the elements (2, 3) and (3, 2) 
in the matrices above. It is then straightforward to see that 
H(k) = (TH(—k))*. Fermionic eigenvectors in the vor- 
tex free sector are thus related to their time-reversed counter- 
parts by conjugation and the reversal of momenta. However 
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(b)Z = 1,J = K = L = 1/^(3) 
FIG. 7: (Color online) Dispersion Relation Eh for (a) the abelian phase, (b) the phase transition (c) the non-abelian phase. 



we again point out that these eigenvectors represent fermionic 
creation and annihilation operators from time reversed sectors 
with different toric code vacuum states. 

The doubled (±E) symmetry of the spectrum means that 
the eigenvalue equation can be written as a cubic polynomial 
equation and analytical expressions for the eigenvalues can be 
obtained, see the appendix. However to calculate the geome- 
try of the phase boundary it is more straightforward to observe 
that within the Abelian phase the minimum gap always occurs 
at k — 0. As this coincides with when A nm = one can cal- 
culate exactly where the phase transition lies by a straightfor- 
ward diagonalization of the £ matrix above with k = 0. The 
eigenvalues are calculated to be 



E = ±2(Z + a^J 2 + K 2 + L 2 ) 



(35) 



where a = —1, or 1 and therefore the minimum energy gap 
is given by \2Z — 2\J 'J 2 + K 2 + L 2 \. The phase transition 
thus occurs at Z = \J J 2 + K 2 + L 2 . This is verified in Fig- 
ure [5] where the minimum energy gap as a function of J = L 
and K was calculated by numerical diagonalization. 

That the minimum energy is obtained at k = also holds 
true along the J = K = L line so long as J = K = L < 

^p-Z. We see therefore that the gap closes and opens lin- 
early as a function of J along this line. The dispersion re- 
lations for the A phase, the critical point , and the gapped 
B-phase along the J = K = L line are shown in Figure [7] 

It is not generally true however that the minimum energy 
occurs when k = 0. This is perhaps most striking when one 
of the parameters J, K, or L are zero, see the appendix and 
Figures [5] and [6] On these planes and outside the radius of 
the sphere we are in a gapless B-phase. In order to open a 
gap we must move off these planes. If for example, K = 
this gapless phase occurs because of two Dirac cones at ±fc. 
Letting K > but with K « J = L opens a mass gap 
in much the same way as the three-body terms do so in the 
original honeycomb model, see Figure|5]and Ref. [TJ In Figure 
[8]we also plot the dispersion relations showing the Dirac cones 
and the opening of the gap with K > 0. Calculations around 
the L — or J = planes reveal similar properties. 



IV. GROUND STATE DEGENERACY ON A TORUS 

The ground state degeneracy on a torus first observed by 
Ref. |2]may be calculated within the Hartree-Fock-Bogoliubov 
formalism in the following way. We have already fixed the 
vacua to be the toric-code states so we know by the Bloch- 
Messiah-Zumino theorem 14 15 that we can write 

v 

|gS>= II n( U <+<W4)l{°' V ' A }>W>- ( 3 6) 

m=l l^m 

where the a^'s are the canonical fermionic raising operators 
gotten from the by performing a singular value decompo- 
sition on the U sub-matrix of the full eigenvector matrix W, 
see for example Ref. 5 The number p gives the number of 
fully occupied un-paired fermions in the ground state and dic- 
tates the fermion number parity. In the extreme case where 
all fermions are unpaired the wavefunction JlL=i a m\ vac ) ' s 
the completely filled Fermi sea. 

On a torus our vacuum state takes the form 



{O, V, A},^ ) ^ 2/) ,{0}) where 1 { Q X> and 1% are the eigen 
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values of the operators and Lq 1 " defined above. For 
the vortex free sector on a torus we can use the Bloch 



(v) 



Hamiltonian ( 3 1 1 to calculate energies, eigenvectors and thus 
the values of ui and u;. In this case the allowed values of k a 
in the various homology sectors on the torus of size (N x , N y ) 
are 9 a + for integer n a — 0, 1, ...N a — 1, where the 

(±1,±1) have values 



four topological sectors, (Zq ,Iq ) 



of 9 a given by 6 a = — )]f L - F° r a l ar g e torus one 
expects that these discretized fc-values becomes so close as 
to give approximately the same ground state energy for each 
homology sector i.e.: 



(37) 



k n—1 



This value comes directly from the assumption that all values 
Ek n associated with the 7^ operators are positive. However, 
cases where some of the 7^ must be associated with nega- 
tive energy solutions do occur and in these cases the ground 
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state energy is raised. This is precisely what happens in the 
( — 1,-1) homology sector. In the non-abelian phase of this 
system any arrangement of the W matrix such that positive 
eigenvalues ( 35 1 at k = are associated with 7^ operators en- 



sures that the matrix U (i.e the upper left quadrant of VV) has 
one column that has just zeros. It is therefore rank deficient 
and, by the Bloch-Messiah-Zumino theorem, the ground state 
has one fully occupied mode (i.e with u = and \v\ — 1) with 
momentum k = 0. 

Assuming we are working on an even-even lattice, the pos- 
sibility of having odd fermion number parity is excluded by 
the fact that we are in the vortex-free sector. This means we 
must switch columns (E7j, V;) <-» (—PV*,PU*) of the eigen- 
vector matrix describing the annihilation and creation of Bo- 
goliubov fermions and perform the singular value decompo- 
sition of the new U matrix. The switching of columns of the 
matrix effectively changes an occupied mode for an empty 
one thereby giving the correct fermionic parity number. How- 
ever, it also raises the ground state energy in this sector to 
E gs = - J2k.n E k,n/2 + Eo,i- If the system is gapped then 
the ground state of this sector is higher than that of the other 
three fully or partially anti-periodic sectors, leading to a re- 
duction in the topological degeneracy. 




-3 -3 



(a)Z = 1, J = L = 1 /%/2 + 0.2, K = 0. 




(b)Z = 1, J = L = 1/V2 + 0.2, K = 0.2 

FIG. 8: (Color online) Dispersion Relation Eh for (a) the gapless 
B-phase with K = 0, (b) the gapped B-phase with K = 0.2. 



V. CONCLUSION AND OVERVIEW 

We have shown that the star-lattice Kitaev model may be 
mapped to a system of fermions hopping on an effective 
Kagome lattice with a Z2 gauge field. The fermionic vacua 
are explicitly shown to be toric code states of the effective 
Kagome spin lattice. The abelian phase of the model is inher- 
ited from the fermionic vacua and time reversal symmetry is 
broken at this level. We have shown that, as in the original 
honeycomb system, there are three distinct phases: a gapped 
Abelian phase, a gapless B-phase and a gapped non- Abelian 
B-phase. The boundary between the phases can be understood 
as a sphere of radius Z in the parameter space of the cou- 
pling strengths J, K and L around the triangles. The gapped 
Abelian phase lies inside the sphere. The gapless B-phase lies 
along the trivial planes J = 0, K = or L = outside 
the sphere. These instances correspond to the cases where 
the underlying lattice is bipartite. Off these planes and out- 
side the sphere the system is in a gapped non-Abelian phase 
with spectral properties similar to those of the gapped non- 
abelian phase of the original honeycomb lattice model. As in 
the original system the ground state degeneracy on a torus can 
be determined by using the Bloch-Messiah theorem. 
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VII. APPENDIX 

In this appendix we will derive the explicit dispersion rela- 
tions for the ground state sectors as well as the exact locations 
of the Dirac cones that occur in the gapless B-phase. The 



(38) 



(39) 



characteristic equation of the 6x6 matrix in ( 32 1 is 

A 6 + aX 4 + bX 2 + c = 0, 

where 

a = -8R(0) - 12Z 2 
b = 16((R(0) + Z 2 ) 2 -2Z 2 (R(k)-Z 2 )) 
c = -64Z 2 (J 4 + K 4 + L 4 + Z 4 +2(T(k) - Z 2 R(k))) 

with R(k) = J 2 cos(k x ) + K 2 cos(k x — k y ) +L 2 cos(fc J/ ) and 
T{k) = K 2 L 2 cos{k x ) + J 2 L 2 cos(k x -k y ) + J 2 K 2 cos(fc y ). 

Equation ( [38] l can be reduced to a degree three polynomial 
by simply substituting X 2 = x. A third degree univariate poly- 
nomial equation c an be solved exactly and the solutions can 
be given as follows^^l. Let q and r be defined as 



q = 



36 -a 2 



9ab - 27c - 2a 3 
54 



(40) 



and the discriminant of this polynomial be p = q 3 + r 2 . Also, 
taking s and t as 



and t = 



- Vp, 



(41) 
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we get three roots xx, X2, X3 as 



x\ 



■''2 



X3 



S + t 
1 



(s + t) - - 
2 y ' 3 



2 V 



. a 
■t)-~ 
' 3 



V3 
2 

2 



(s - t)i 
(s - t)i. 



(42) 



Now, for p > 0, two of these are complex conjugate roots. 
However, as all the eigenvalues are known to be real we may 
disregard this case. Forp < 0, we let 9 = arccos( J- — ), and 



the three real roots can be given as 



X\ 
Xl 

x 3 



2sf=q cos(-) - - 

9 + 27T 

2V-? cos( — - — ) 



= 2 A 



- ,6> + 47T. 

-q cos( — - — ) 



(43) 



The real roots of a generic cubic univariate polynomial equa- 
tion have recently been classified in term s of its coefficients 
using real algebraic geometry techniques'^. For p = 0, 
r = ±y — g 3 and so X2 = £3. For all the cases we have 
observed the dispersion relation (the minimum energy band) 
is given by E k — A 2 = y^. 

Gapless phases in the system may be identified by setting 
c = or equivalently setting det(-ff) = 0. From this we 
obtain a quadratic equation in the variable z = Z 2 : 

z 2 + 2R(k)z + 2T(k) + J 4 + K 4 + L 4 = 0. (44) 

with generic solutions 



R(k) ± ^fR{k) 2 - [2T(k) + J 4 + K 4 + L% 



(45) 



If k = we have R{0) 2 = [2T(0) + J 4 + K 4 + i 4 ] and 
the quantum phase transition can be seen to occur on the sur- 
face of the 3 -sphere Z = y/R(0) in the space of the coupling 
parameters J,K and L, see Figure [6] 

Other interesting cases occur when either of the parameters 
J,K or L vanish. In these situations we have 



i?(fc) 2 



K 4 



hi 1 



L 4 

-2 



L sm(ky) 



- [2T(k) + J 4 

= i\K 2 sin(fc a; - 

= i[J 2 sin(fc x ) + L 2 sin(fc 2/ )], 

= i[J 2 sm.(k x ) + K 2 sm{k x — k y )], 



(46) 



J = 
K = 
L = 



which, because the parameter Z is real, must vanish. This can 
be used to set conditions on parameters k x and k y . In the cases 
above if the two non-zero terms are set to be equal (for exam- 
ple L — K = M, J = 0), the three equations above reduce 
to (i) k x = 2k y (ii) k x = — k y (iii) 2k x — k y respectively. 
Substitution into the remaining part Z = \J R(k) gives the 
zero energy solutions at 



(k x 1 k y ) 



J = 0,L = K = M 
K = 0, J = L = M 
L = 0, J = K = M 



±(21,1) 
±(1,-1) 
±(1,21) 

where I = cos" 1 (Z 2 /2M 2 ). If M < l/y/2 we see that I is 
complex but as we have already shown these values of J, K 
and L corresponds to the gapped phase inside the sphere Z = 
y/R(0). The zero energy solutions described here form the 
sharp points of the Dirac cones at these values of k x and k y . A 
gap may be opened in the spectrum at these conical points by 
letting each of the parameters J, K and L have non-vanishing 
values such that the condition \J R(0) > Z is fulfilled. 
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